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ABSTRACT 

We perform numerical simulations of solid particle motion in a shearing box model 
of a protoplanetary disc. The accretion flow is turbulent due to the action of the 
magnetorotational instability. Aerodynamic drag on the particles is modelled using 
the Epstein law with the gas velocity interpolated to the particle position. The effect 
of the magnetohydrodynamic turbulence on particle velocity dispersions is quantified 
for solids of different stopping times ts, or equivalently, different sizes. The anisotropy 
of the turbulence is reflected upon the dispersions of the particle velocity components, 
with the radial component larger than both the azimuthal and vertical components for 
particles larger than ~ 10 cm (assuming minimum-mass solar nebula conditions at 5 
AU). The dispersion of the particle velocity magnitude, as well as that of the radial and 
azimuthal components, as functions of stopping time, agree with previous analytical 
results for isotropic turbulence. The relative speed between pairs of particles with the 
same value of ts decays faster with decreasing separation than in the case of solids 
with different stopping time. Correlations in the particle number density introduce a 
non-uniform spatial distribution of solids in the 10 to 100 cm size range. Any clump 
of particles is disrupted by the turbulence in less than one tenth of an orbital period, 
and the maximally concentrated clumps are stable against self-gravitational collapse. 

Key words: accretion, accretion discs - Solar system: formation - planetary systems: 
protoplanetary discs - MHD - turbulence 



1 INTRODUCTION 

In order for planetesimals to form in a circumstellar disc, the 
parent dust grains must meet a variety of kinematic condi- 
tions if they are to grow in mass and size. The formation of 
kilometre-sized solid bodies by gravitational instability re- 
quires, through the Toomre criterion, that the dust spatial 
density be pd > /nG, where Q is the disc angular veloc- 
ity and G the gravitational constant. Such densities could be 
achieved in a mid-plane dust layer after sedimentation has 
occurred. On the other hand, growth by coUisional sticking 
may be accomplished if the grain relative speeds acquire cer- 
tain values: experiments have shown that, for micron-sized 
grains, collision speeds less than ~0.2 m/s lead to aggregates 
with a fractal geometry (Blum & Wurm 2000), and above ~ 
1 m/s disruption of the colliding agglomerates takes place. 
While binary collisions of the smallest dust grains are 
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caused by Brownian motion, other sources of relative ve- 
locities in a protoplanetary nebula may come into play as 
the interaction between solids and gas decreases. Differen- 
tial sedimentation, radial drift and turbulence can provide 
the necessary kinetic energies to drive the coUisional motion 
of solids, in addition to changing their spatial distribution 
throughout the nebular disc. In this paper we focus on the 
effect of a turbulent gas flow on solid particle velocities and 
spatial arrangement. 

The formalism of Volk et al. (1980) makes it possible 
to calculate the velocity dispersions of particles in isotropic 
turbulence. The effect of the fluctuating gas velocity on indi- 
vidual solids is approximated by a sum over two types of tur- 
bulent eddies: those with spatial frequencies k and turnover 
times tk > ts, where ts is the particle stopping time, and 
which merely advect the particle; and those with tk < ts, 
which decay before the particle has had time to cross them. 
The result is a velocity dispersion that decreases with parti- 
cle stopping time as (f2ts)~^^^, with a sharp steepening for 
particles with fits ~ a few times 10~^. Under minimum mass 
solar nebula conditions at 5 AU, this corresponds to parti- 
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cle radii of tens of centimetres. The method is nicely sum- 
marised and applied by Cuzzi & Hogan (2003) to chondrule- 
sized objects. 

The accumulation of a large number of solids in small 
regions of the turbulent protoplanetary disc can lead to feed- 
back effects on the gas flow if the solid mass density is greater 
than the local gas density. Such particle agglomerations can 
also lead to dust growth if the associated collision rates 
are conducive to sticking. One recently suggested clumping 
mechanism is a streaming instability that arises due to the 
relative motion between solids and gas in a Keplerian disk 
(Youdin & Goodman 2005). The feedback of solids onto the 
gas can generate exponential growth of particle density per- 
turbations as well as turbulence, and particle clumps can 
develop overdensities up to a factor of ~ 50 relative to the 
mean gas density (Youdin & Johansen 2007, Johansen & 
Youdin 2007). 

The process of turbulent concentration (Maxey 1987) 
has been explored by Cuzzi et al. (2001), who use the ratio of 
local particle density to its global average as the measure of 
a multifractal distribution of concentrated chondrules. This 
distribution resembles that of the dissipation of the turbu- 
lent kinetic energy on the Kolmogorov scale (Chhabra 1989). 
The connection between turbulent concentration and tur- 
bulent dissipation suggests that the local structure of the 
chondrule concentration field may be independent of flow 
properties. For larger particles, an accurate characterisation 
of their number statistics in the turbulent circumstellar en- 
vironment is still lacking. 

The recognition that turbulence due to the magnetoro- 
tational instability (MRI; Balbus & Hawley 1991; Hawley, 
Gammie & Balbus 1995, henceforth HGB) can provide the 
viscous stresses necessary to allow inward accretion in discs, 
has prompted investigations of different dynamical aspects 
of dust evolution in the presence of magnetohydrodynamic 
(MIfD) turbulence. Planetesimal formation by gravitational 
instability was addressed by Johansen et al. (2006, hence- 
forth JKH), whose numerical studies concluded that groups 
of at least 1000 particles in their simulations could become 
gravitationally bound, assuming specific values of disc col- 
umn densities (150 and 900 g cm~^) and a ratio of disc scale 
height to radius of 0.04. Fromang & Nelson (2005) examine 
the accumulation of individual solids due to vortical struc- 
tures in a global protoplanetary disc model. They follow the 
radial migration of up to 3000 bodies, whose drag interac- 
tion with the gas corresponds to objects of approximately 1 
m in size. Those particles able to migrate inwards do so at 
different rates, with reductions in semi-major axes ranging 
from a factor of 1.6 in approximately 50 orbits to a factor 
of 2.2 in 200 orbits. The particles that become trapped in 
vortices maintain an approximate constant radial position 
until the end of the simulation. 

In the following we present a study of the kinematics 
of solid particles in a local MHD model of a protoplanetary 
disc. In Section 2 we describe the numerical method em- 
ployed to model the accretion flow and the particles. Section 
3 contains results obtained from measurements of particle 
velocities and spatial distribution. In Section 4 we discuss 
their significance in the context of particle growth in a pro- 
toplanetary nebula, and concluding remarks are presented 
in Section 5. 



2 METHOD 

2.1 Numerical Method 

We use a three-dimensional version of the ZEUS code (Stone 
& Norman 1992a;b) in which the shearing box model of 
HGB is implemented. To simulate the solid particles, we 
have added a module that solves their equation of motion, 

1 o 

f= (vp - Vg) - 2n X Vg -I- 3f7 XX (1) 

ts 

where f is the force per unit mass on the particle, the first 
term on the right-hand side is the Epstein drag law for spher- 
ical particles, ts is the particle stopping time, Vp is the parti- 
cle velocity, Vg is the velocity of the background fiow at the 
particle's position, and f2 = (0, 0, fi) is the angular velocity 
of the box. The last two terms of Eq. ^ are the Coriolis 
and tidal forces, respectively, that arise as a result of the 
non-inertial character of the shearing box {x is the distance 
in the radial direction x). The particle module takes as in- 
put the gas velocities v calculated by ZEUS on the faces of 
the computational grid cells. These velocities are obtained 
by solving the ideal MHD equations for the shearing box, 
together with an isothermal equation of state, 

^ + V • (pgv) = (2) 

^+v.Vv = -1 V f P + ^1 + (^-^)^-2»xv+3»-xx 

(3) 

^ = V X (v X B) (4) 



P = Pgc? (5) 

where Pg is the gas density, B is the magnetic field, P is 
the gas pressure, and Cs is the isothermal sound speed. The 
module interpolates each velocity component to the posi- 
tion of each particle inside the respective computational 
grid cell, using a trilinear interpolation algorithm (Press et 
al. 1992). Equation Q is solved via a second-order Runge- 
Kutta method. The back reaction of the particles on the 
gas, inter-particle interactions, and particle coupling to the 
magnetic field are not included. 

Despite the likely presence of dead zones in protoplane- 
tary discs as a result of non-ideal MHD processes (Gammie 
1996), in this first study we neglect the effect of resistivity 
and assume that our local disc model represents an active 
region in the surface layers or at a large orbital radius (Sano 
et al. 2000). 

We adopt the values used in HGB for the initial density, 
initial pressure, angular velocity, plasma parameter and box 
dimensions: po = 1, Po = 10"*^, Cs = = 10"^, /? = 400 and 
H X 2nH X H (where H is the disc scale height) , respectively. 
The initial magnetic field is vertical but, unlike HGB, it 
varies sinusoidally in the radial direction, in such a way that 
the net magnetic fiux is zero. The grid resolution used is 
84 X 180 X 84 zones. 
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Simple tests of the particle integrator involved the evo- 
lution of a single solid particle in a zero- velocity background 
flow, with various values of the stopping time ts which, in 
the Epstein regime, is given by (Weidenschilling 1977). 



CsPg 

where a is the particle radius and pp the particle solid den- 
sity. Given an initial position and velocity, the response time 
of the particle to the drag force was found to adjust ac- 
curately to the expected analytical behaviour for a range 
of values of the dimensionless stopping time Qta between 
2x10"-^ and 100. 

2.2 Simulations 

We performed two sets of simulations in which the shearing 
box is first evolved for 10 orbits (1 orbit=27r/f2) to ensure 
that turbulence due to the MRI has developed. At this time, 
A'^ particles are introduced in the flow with random initial 
positions and with initial velocities set to zero. In simulation 
set A, A'^ = 80, 000, and the total number of particles is 
divided in eight groups of 10,000 members. Each group is 
given a specific value of The values used were 100, 30, 
10, 5, 2, 0.2, 0.02 and 0.002. The system was evolved for a 
further 100 orbits. 

Simulation set B, which was used to study the par- 
ticle spatial distribution, consists of four groups of A = 
4, 704, 000 particles, with each group having one of ^lts=2, 
0.2, 0.125 and 0.05. The code was run for a further 30 orbits. 



Also shown as plus signs are the ratios Sv = cFv/'^v of the 
velocity magnitude dispersions a^'^ = •\J {{vp^g — {vp,g))^) 
(notice that 5*^5 ^ + Sy + Sf). The dotted curve corre- 
sponds to the analytical solution of Volk et al. (1980) and 
Cuzzi et al. (1993) for Sv The solid and dashed curves cor- 
respond, respectively, to the radial and azimuthal solutions 
of Youdin and Lithwick (2007). With no gas drag, the par- 
ticles would follow epicyclic paths with the ratio a^^/a^^ 
of the radial and azimuthal velocity dispersions equal to 2. 
Larger ratios a^^/u^^ of around 2.3 measured in the MHD 
calculation are due to weak drag forces from the turbulence. 
Separate numerical integration of the particle equation of 
motion ([1} with a sinusoidally time-varying gas velocity re- 
produces the anisotropy, given variation periods of one to 
two orbits. Fourier analysis of the gas velocities at the parti- 
cle locations in the MHD calculation indicates the strongest 
modes do have periods longer than one orbit. Our neglect 
of the vertical component of gravity appears to have little 
effect on the ratios between particle velocity components, 
as similar anisotropy was observed in stratified shearing box 
calculations (Carballido et al. 2006), in which for the f2ts=10 
case, (jJJ^ 2.5ctJ^ . 

It is important to remark that the analytical models of 
Fig. 2 were originally derived assuming an isotropic turbu- 
lent velocity field, so in principle one should not expect that 
they be an accurate description of particle motion in MRI 
turbulence. Nevertheless, the particle data obtained from 
the simulations conform reasonably well to the models. This 
will be discussed below. 



3 RESULTS 

3.1 Viscous evolution of disc 

Figure 1 shows the evolution of the magnetic {dotted line) 
and hydrodynamic {dashed line) volume-averaged stresses 
associated with our disc model, during the longest run per- 
formed. The stresses are normalized to the initial pressure. 
Their sum {solid line) gives an estimate of the volume- 
averaged a parameter of classical disc theory. The average 
value of a is 8 X 10~^, which falls in the range 10~* - 10"^ 
inferred for T Tauri stars (e.g. Hartmann et al. 1998). 

3.2 Velocity dispersions 

Simulation set A has been used to calculate the time- 
averaged distributions of the particle velocity components. 
The resulting dispersions, crjj. (i = x, y, z), normalised by the 
corresponding gas velocity dispersions erf., have been mea- 
sured as functions of particle stopping time. To calculate 
this ratio, Si = a^. /erf. , the gas velocities were computed at 
the position of the particles, at each snapshot and for each 
value of Qts, and the effect of shear was subtracted from the 
particle and gas y- velocities. The ratios Si thus obtained 
were then averaged over the last 10 orbits of the simulation, 
a period in which the largest particles (those with f2ts = 100) 
have reached dynamical equilibrium with the gas. 

The results are summarised in Fig. 2, where the solid- 
to-gas ratios of the x {squares), y {triangles) and z {filled 
circles) components of the velocity dispersions are plotted. 



3.3 Pair-wise relative velocities 

Fig. 3 shows the magnitude Vj-^i of the relative velocity of 
pairs of particles that are located inside a cubic sub- volume 
in the box, in simulation set A, as a function of inter-particle 
separation r (in units of the box size H). The velocity has 
been normalised to the gas sound speed. The sub-volume 
has dimensions 1x1x1, and extends along the y direc- 
tion between the values y = 2.64 and y — 3.64. At each 
snapshot, a sample of 700 particles, all with the same stop- 
ping time, is randomly selected in the sub-volume, and the 
relative speed for each of the 244,650 possible pairs in the 
sample is calculated. The relative separation between pair 
members is binned, with the bin size equal to the a;-length 
of a grid cell. The values of Vrei are then averaged in each 
bin over a duration of 7 orbits. The data shown correspond 
to the stopping times f2ts = 100 {triangles), 10 {squares), 0.2 
{diamonds), and 2xl0~^ {filled circles). To calculate Uroi, 
the contribution of the shear has been subtracted from the 
y component of the particle velocities. 

It can be seen that at small separations the relative 
speed decreases appreciably. This is to be expected, since 
grains that are equally coupled to the gas move coherently 
with small turbulent eddies. In the limit of particle stop- 
ping times less than the lifetime of the smallest eddies, the 
laminar flow on scales below the dissipation scale does not 
contribute to the random relative grain velocity. Formally, 
this can be expressed for grains of stopping times ts-^ and 
as (Weidenschilling 1984) 
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Figure 1. Time evolution of the stresses in the local disc model. The stresses are normalized to the initial gas pressure. The sum of the 
magnetic and hydrodynamic stresses (solid line) provides an estimate of the viscous a parameter, whose average value is 8 X 10"'^. 



where u and r are the velocity and lifetime of the smallest 
eddy, respectively. 

A similar calculation as that of Fig. 3 is presented in Fig. 
4, but for pairs of particles with different stopping times. 
Each symbol represents a different pair of stopping times 
(the error bars have been ommited for clarity, but give an 
uncertainty of ~ 10% for the data represented by the plus 
signs, the diamonds and the triangles, and of ~ 5% for the 
filled circles). The relative speeds do not decrease as steeply 
with decreasing separation as in the case of like particles. 
Again, from the order-of-magnitude estimate ([7} , Vrai should 
acquire a finite value in the limit of small stopping times. For 
the particles considered, this value is a significant fraction 
of the sound speed (~ 12 — 15%). 

It is worthwhile to point out that as the inter-particle 
separation decreases, the effects of numerical resolution may 
become important for the purpose of calculating the rela- 
tive velocities. The size I of the smallest turbulent eddies, 
to which the smallest particles are coupled, varies with vis- 
cosity 11 as £ ^ z/"^''*. The viscosity in turn can be affected 
by the grid resolution (HGB) . 



3.4 Spatial distribution 

In order to improve the statistics of particle counts, we have 
taken as the unit of volume a group of 3 x 3 x 3 grid cells, 
referred to as a 3-macrocell. In this way, the average oc- 
cupation number is 100 particles per 3-macrocell, with an 
associated standard deviation of 10%. 

The placement of the dust particles in the shearing box, 
at time i=10 orbits, is random, and therefore follows a Pois- 
son distribution. This may no longer be true at subsequent 
times, when the turbulence has mixed the particles in such 
a way that the initial stochastic distribution has been lost. 
This situation is apparent in Fig. 5, where the left panel 
shows the cumulative distribution function, as a solid line, 
of the particle number n per 3-macrocell at the moment 
when the particles are introduced in simulation set B, for 
the case. The dotted line is the corresponding cumu- 

lative Poisson distribution with the same mean ((n) = 100) 
as the n distribution. The two are alike, indicating that the 
initial particle positions are indeed drawn from a random 
sample. However, the right panel, which corresponds to t—30 
orbits, makes it evident that the spatial distribution of the 
particles deviates from a Poisson distribution. The vertical 
dashed line in each panel indicates the point at which the 
distance d between the two distributions is greatest. This 
distance constitutes a measure of the deviation of one data 
set from the other. In the particular case shown in Fig. 5, 
at t=10 orbits the maximum distance is d=0.026, whereas 
at t—30 orbits, d—0.55. The situation is similar at all other 
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Figure 2. Ratio of particle to gas velocity dispersion, as a function of particle stopping time. Data shown correspond to the radial 
(squares), aaimuthal (triangles) and vertical (circles) components of velocity, as well as to the velocity magnitude (plus signs). The 
curves represent analytical models derived for isotropic turbulence. 



times t >0, and for the three other hydrodynamic couphngs 
of set B. 

To measure the amount of solid particle clumping in the 

shearing box, we calculate the number of particles in each 
macrocell and divide by its mean: 



Fig. 6 shows distributions of C for each particle-gas 

coupling, at t=20 orbits. These distributions are typical 
throughout the duration of the simulation. Also shown for 
comparison is the distribution of gas clumping, Pg/(pg). No- 
tice that the solids exhibit more over-densities with respect 
to the mean ((C) — 1) than under-densities, showing a con- 
siderably high rmmber of rnacrocells that have values of C 
above the average occupation rmmber (n). The highest C 
attained is Cmax=26.5, by the particles with nts=0.125, at 
t =27.1 orbits. A history of Cmax is plotted in Fig. 7, where 
the time axis has been shifted to the instant at which the 
particles are introduced. 

The time-averaged C distributions are shown in Fig. 8. 
Note the difference in the shapes of the distributions com- 
pared to those in Fig. 6, particularly for Qta = 0.2, 0.125 
and 0.05 at C = 0. This can be understood by consider- 
ing that the time-averaging procedure effectively assigns a 
non-zero number of particles to each 3-macrocell, which at 
a specific snapshot (but not all) is likely to be empty. The 



non-Poissonian character of the particle spatial distribution 
is also evident in Fig. 8, as the distributions possess ap- 
preciable "tails" extending in the direction of increasing C. 
This is normally associated with spatial correlations in a 
system of particles in which their positions are not statis- 
tically independent (Landau & Lifshitz 1980). In the case 
of particles immersed in turbulence, these correlations arise 
from the intermittency of turbulent mixing. The correlations 
arc manifested by an increase in the variance of the particle 
counts, with respect to a purely random distribution. 

Positive correlations can induce clustering at some 
range of spatial scales. One way to quantify this clustering 
is by means of the pair correlation function 



/A N (n(r)n(r + Ar)) 
^(^'^^ = (n(.))^ 

where the particle number density n is evaluated at two 
points separated by the distance Ar. The measurement of 
this quantity is common in the study of cloud droplet clus- 
tering in the presence of turbulence in the Earth's atmo- 
sphere (Kostinski & Jameson 2000). In order to calculate r] 
for the population of solids in our shearing box, we proceed 
as follows. Given the anisotropic character of the turbulent 
eddies present in the box (HGB), 77 is measured indepen- 
dently along each spatial direction i (=1,2,3, or x,y,z): 
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Figure 3. Relative speed (in units of the sound speed) as a function of inter-particle separation (in units of scaleheight), for pairs of 
particles with the stopping times shown on the upper left-hand corner. The number of particles per stopping time is 700. Data have been 
binned with respect to r, with a bin size equal to the size of a grid cell in the x direction (Ri O.OliJ). The values of ii^ei decrease with 
small values of the separation. 



{n{xi)n{xi + Axi)) 



(10) 



Furthermore, to obtain a sufficient number of points 
along each direction, at which to evaluate n, we consider 
2-macrocells instead of the 3-macrocells used to calculate 
the clumping factor C. To avoid the effect of the periodic 
boundary conditions in the y and z directions, as well as 
the quasi-periodic conditions along the radial direction x, 
the separation Axi between 2-macrocells is constrained to 
be, at most, half of the corresponding box length. Thus, for 
example, to evaluate rjx at separations of one macrocell, a 
horizontal row of grid cells is chosen, say one with end cells at 
x = —0.5 and x = 0.5, and with fixed y and z coordinates. 
The particle number density n is measured at the first 2- 
macrocell in this row, and also at the next macrocell. The 
product of these two values of n is stored, and the process 
is repeated starting at the second macrocell. If m is the 
number of 2-macrocells in the row, the process ends with 
the product of n evaluated at macrocells m/2 — 1 and m/2. 
The average of all the products is calculated and inserted in 
the right-hand side of Eq. ((TD}. 

The row of grid cells used in the calculation just de- 
scribed is located at specific y and z, and it is only one 
of 84 X 180 = 15, 120 similar horizontal rows. With this 
number of values for rix{Ax), an average over all rows at 



all y and z gives the volume- averaged pair correlation func- 
tion {rix{Ax)). The functions {r]y[Ay)) and {r]z[Az)) are ob- 
tained in an analogous manner. 

The volume-averaged pair correlation functions for the 
f2ts=0.05 particles, at the time of their introduction in the 
flow {t = 10.0 orbits), are identically zero, reflecting the sta- 
tistical independence of the initial particle numbers. This 
also holds for the cases nfs=0.125, 0.2 and 2. The mea- 
sured {r]i) at the end of the simulation are shown in Fig. 9. 
The columns correspond to the x, y and z directions, and 
the rows represent, from top to bottom, the stopping times 
^Its —2, 0.2, 0.125 and 0.05. The error bars correspond to 
the data dispersion at each separation Axi. 

One feature is immediately noticeable from Fig. 9: the 
vertical correlation function (third column) is very nearly 
flat for the four values of fits. This would seem to indicate 
that particles remain almost evenly mixed in the vertical 
direction, without deviating much from the initial random 
distribution. The particles show concentrations over radial 
distances of approximately O.IH (for Qts—2) and 0.05H (for 
the other stopping times). In the azimuthal direction the 
respective correlation lengths are ~ H and ~ 0.3H. This 
range of distances represents between 5% and 16% of the 
box size. 

The largest correlation lengths belong to the Qta=2 
particles, which also show the lowest average concentration 
(Figs. 7 and 8). Even taking the highest value of Cmax for 
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Figure 4. Same as Fig. 3, but for pairs of particles with different stopping times. 
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Figure 5. Cumulative distributions of the number of particles per 3-macrocell, for Cits = 2. The left panel represents the initial 
distribution at t=10 orbits. The solid line corresponds to the particle distribution n, while the dotted line is the corresponding Poisson 
distribution with the same mean. The vertical dashed line indicates the value of n at which the vertical separation d between the two 
distributions is greatest. The right panel shows that the particle distribution is not Poissonian at t=30 orbits. This is true in general for 
t >0 and the three other values of fits ■ 



this stopping time (~21), the surface density reached in a 
clump of dimensions (rix) x {rfy) is w 2.1 x 10*i?~^, com- 
pared to 2.65 X W^H~^ for the most highly concentrated 
particles, those with fits =0.125. 

In order to constrain the lifetime of typical particle as- 
sociations from the simulation data, it is necessary to refine 
the working definition of a clump. Here a clump is under- 



stood to be a group of particles that has a size smaller than 
the length 6 of the diagonal of a 3-macrocell. The size of 
the group is determined by calculating a clump "diameter" , 
based on the maximum and minimum values of the particle 
coordinates. 

For each value Qts=0.2, 0.125 and 0.05, the clump with 
the highest particle density was followed during subsequent 
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Figure 6. Distribution of the clumping number C = n/{n) at t =20 orbits, for all particle stopping times. All solids exhibit an 
overabundance of 3-macrocells for which n is larger than the mean. All stopping times are more clumped than the gas (dash-tripple 
dotted line). 



snapshots. In all cases, the particles become detached from 
the original clump after At = 0.1 orbits, in such a way 
that a calculation of the new diameter yields a value larger 
than 5. The fact that At is the interval between snapshots 
indicates that the lifetimes are less than what is possible 
to probe using the available data. The situation is similar 
if by clumps it is meant groups with horizontal dimensions 
corresponding to the particle correlation length. 



4 DISCUSSION 

The measurement of solid particle velocity dispersions in 
MRI turbulence allows to test the analytic model of Volk 
et al. (1980) and Cuzzi et al. (1993), which adopted a Kol- 
mogorov energy spectrum. The turbulence that develops in 
the shearing box setup that we use follows roughly a Kol- 
mogorov scaling (HGB), and the dispersion of the particle 
velocity magnitude is generally consistent with the previous 
results. 

The radial and azimuthal velocity dispersions were 
found to conform to the analytical results of Youdin and 
Lithwick (2007), who use a turbulence model in which gas 
velocity correlations {v^Vy) are weak compared to {v'^ +Vy)- 
This is indeed the case in the MRI-generated turbulent fiow 
present in our shearing box. 

Recent analytical calculations have provided closed- 
form expressions for the relative velocity between solids of 
different sizes, as a function of stopping time, in isotropic 



turbulence (Ormel & Cuzzi 2007). Those results validate the 
use of Eq. ((Tjl in a regime where the stopping times of both 
particles are smaller than the lifetime t of the smallest ed- 
dies. From a measurement of the correlation time of the gas 
velocities, we find t ~ 0.1 orbits. Therefore, Eq. ((Tjl will not 
apply to all the stopping times shown in Fig. 4, since the 
stopping times for the S74 = 10 and f2ts=100 particles are, 
respectively, ~1.6 and ~16 orbits. Nevertheless, a separate 
calculation shows that relative velocities between pairs of 
the three smallest stopping times (Slfs=0.002, 0.02 and 0.2, 
which correspond to ~ 3x 10""', 3x 10~^ and 3x 10~^ orbits) 
are well described by Eq. (O at short separations, with the 
lowest velocity {vj-^i ~ 0.09cs) corresponding to pairs with 
fits =0.002 and 0.02. 

For ts > tL, with tL being the overturn time of the 
largest eddies (typically taken as one orbital period), the 
model of Ormel & Cuzzi (2007) gives 



2 _ 2 



1 + nts, 1 



(11) 



where Vgas is the velocity of the turbulent flow, and it is 
assumed that ta^ > ts2 whitout loss of generality. Direct 
measurement from our simulation data yields Vg^s ~ 0.12cs, 
and so for the two largest stopping times studied (S7ts=10 
and r2ta=100) we obtain Vj-d ~ 0.04cs, in close agreement 
with the numerical results of Fig. 4 (filled circles). 

The two-point correlation function calculated in Sec- 
tion 3.4 gives an estimate of the length scale over which 
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Figure 7. Time history of the maximum clumping factor Cmax for the particle stopping times used in run B. The time marked as t = 
orbits corresponds to the instant at which particles are introduced in the flow. 



particle clustering occurs. In a MMSN scenario at 5 AU, 
in which the scale height H of the disc is ~ 4 x 10^^ cm, 
the minimum correlation length of the particle number den- 
sity would be ^ ~ 0.05H ~ 2 x 10^ km. A recent model of 
chondrule formation, in which molten chondrules come to 
equilibrium with the gas evaporated from other chondrules, 
indicates that these particles must have originated in regions 
larger than 6,000 km in radius, or ^ ~ 6 x 10~^H (Cuzzi & 
Alexander 2006). The smallest particles in our study of the 
spatial distribution correspond roughly to solids with a ra- 
dius of 20 cm, two orders of magnitude larger than typical 
chondrule sizes. It would be necessary to perform additional 
calculations with particle stopping times that correspond to 
this size range, in order to further constrain the value of ^. 
However, it is unlikely that current computational resources 
would allow for the necessary resolution to probe such small 
correlation lengths. 

We have estimated the Jeans radius Rj of the maxi- 
mally concentrated clump for each stopping time, following 
calculations by JKH. Using similar disc parameters (column 
density E = 900 g cm~^, average dust-to-gas density ratio 
= 0.01, and scale height-to-radius ratio H/r = 0.04) and 
a turbulent diffusion coeflicient (5t = a = 8 x 10~^ inferred 
from our shearing box model, the values of Rj are smaller 
than the typical clump size by factors between ~9 and '^250. 
The clump radius that is closest to Rj corresponds to the 
nta=2 particles. Under these assumptions, our dense par- 
ticle associations would not be subject to a gravitational 
collapse. 



Even though the stopping times used in this calculation 
and those studied by JKH are different, we can still compare 
our results for flts=2 with their Qta=l case. The main dif- 
ference resides in the particle velocity dispersion, which we 
measure as 0.036ca , approximately 50% larger than they ob- 
tain. If our maximally concentrated Qts=2 group were to be 
gravitationally unstable, the velocity dispersion would need 
to be a factor ^ 2.5 smaller. We also note that the radial 
drift incorporated by JKH leads to particle concentrations 
that are stronger than in the cases presented here (Johansen 
et al. 2007), making gravitational collapse more likely. 



5 CONCLUSIONS 

We have performed numerical MHD calculations of solid 

particle velocities and spatial distribution in a turbulent pro- 
toplanetary disc, in the context of the shearing box model. 
The results suggest that, even though the turbulent flow 
has an anisotropic spectrum, the dispersion of the particle 
velocity magnitude (normalised by the corresponding dis- 
persion for the gas) as a function of particle stopping time 
is consistent with analytical models that assume isotropy. 
Nevertheless, the dispersions of the individual velocity com- 
ponents for each stopping time are different. In the case of 
particles with fits }^ 1, the radial component of the solid- 
to-gas velocity dispersion ratio is larger than the azimuthal 
and vertical components typically by factors of 2.3— ~ 
2.5. The relative speed between pairs of particles with the 
same stopping time decreases towards small values at small 
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Figure 8. Time-avcraged distribution of the clumping number C, for all particle stopping times. The tails of the distributions, although 
not particularly long, are indicative of the non-Poissonian nature of the number density counts. 



separations, whereas for pairs of different stopping times the 
relative speed reaches a finite value that depends on the sizes 
of the pair members. 

The clustering induced on the particle spatial distribu- 
tion by the turbulent flow gives rise to short-lived clumps. 
Under conditions of a minimum-mass solar nebula, these 
clumps are at least nine times larger than their respective 
Jeans radii, as a result of high particle velocity dispersions. 
Gravitational collapse would therefore seem to be an un- 
likely mechanism in the formation of more massive objects, 
under the conditions that we assumed. 

Future calculations should increase the numerical res- 
olution used in the measurement of the relative velocities 
between particles, in order to improve fittings to existing 
analytical models. 
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Figure 9. Pair correlation function for all particle stopping times, at the end of simulation set B, The left, centre and right columns 
correspond to the x, y and z directions, respectively. First row: Qts = 2; second row: Htg = 0,2; third row: Qts = 0,125; fourth row: 
Qts = 0.05, The vertical correlation functions are very nearly zero, indicating that particles are evenly mixed in that direction. 
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